Demonstration for Comparing the MSIS model versions in GEODYN¶
This notebook shows the output for running GEODYN with the three latest versions of the MSIS series of upper atmospheric models: MSIS_86, MSIS_00, and MSIS_2.0.
The challenges that were overcome to include these three models into versions of GEODYN are detailed [add page detailing the challenges]
Input the Parameters of the Starlette-SLR GEODYN run:¶
input the parameters that make up the starlette-SLR run type in GEODYN
further specify the runs to the different density model output.
[1]:
%load_ext autoreload
%autoreload 2
##################################
# INPUT PARAMETERS:
##################################
sat = 'st'
arc = '030914_2wk'
grav_id ='goco05s'
local_path = '/data/analysis/starlette_analysis/'
Accel_Status = 'acceloff'
SAT_ID = 7501001
###################################
#
# CHOOSE Models to compare: #
m1 = 'msis86' #
m2 = 'msis00' #
m3 = 'msis2' #
###################################
##################################
# PATH TO DENSITY MODEL RUNS:
##################################
msis86_model = 'msis86'
path_to_msis86 = '/data/runs_geodyn/'+sat+'/results/'+msis86_model+'/'+ msis86_model+'_'+ Accel_Status +'/'
msis2_model = 'msis2'
path_to_msis2 = '/data/runs_geodyn/'+sat+'/results/'+msis2_model+'/'+ msis2_model+'_'+ Accel_Status +'/'
msis00_model = 'msis00'
path_to_msis00 = '/data/runs_geodyn/'+sat+'/results/'+msis00_model+'/'+ msis00_model+'_'+ Accel_Status +'/'
[2]:
if m1 == 'msis86':
path_to_m1 = path_to_msis86
elif m1 == 'msis00':
path_to_m1 = path_to_msis00
elif m1 == 'msis2':
path_to_m1 = path_to_msis2
if m2 == 'msis86':
path_to_m2 = path_to_msis86
elif m2 == 'msis00':
path_to_m2 = path_to_msis00
elif m2 == 'msis2':
path_to_m2 = path_to_msis2
if m3 == 'msis86':
path_to_m3 = path_to_msis86
elif m3 == 'msis00':
path_to_m3 = path_to_msis00
elif m3 == 'msis2':
path_to_m3 = path_to_msis2
[3]:
import sys
sys.path.insert(0, '/data/analysis/util_funcs/py_read_geodyn_output/')
from a_ReadGEODYN import Read_GEODYN_func
Load the simulation output for each density model¶
[4]:
(AdjustedParams_m1,
Trajectory_m1,
Density_m1,
Resids_m1) = Read_GEODYN_func(arc,
sat,
SAT_ID,
grav_id,
local_path,
path_to_m1,
False,
2003,
'SLR')
The base file name for this arc is: st030914_2wk.goco05s
File exists: iieout
/data/runs_geodyn/st/results/msis86/msis86_acceloff/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
/data/runs_geodyn/st/results/msis86/msis86_acceloff/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
/data/runs_geodyn/st/results/msis86/msis86_acceloff/DENSITY/st030914_2wk.goco05s
Loading data...
Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded
[5]:
(AdjustedParams_m2,
Trajectory_m2,
Density_m2,
Resids_m2) = Read_GEODYN_func(arc,
sat,
SAT_ID,
grav_id,
local_path,
path_to_m2,
False,
2003,
'SLR')
The base file name for this arc is: st030914_2wk.goco05s
File exists: iieout
/data/runs_geodyn/st/results/msis00/msis00_acceloff/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
/data/runs_geodyn/st/results/msis00/msis00_acceloff/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
/data/runs_geodyn/st/results/msis00/msis00_acceloff/DENSITY/st030914_2wk.goco05s
Loading data...
Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded
[6]:
(AdjustedParams_m3,
Trajectory_m3,
Density_m3,
Resids_m3) = Read_GEODYN_func(arc,
sat,
SAT_ID,
grav_id,
local_path,
path_to_m3,
False,
2003,
'SLR')
The base file name for this arc is: st030914_2wk.goco05s
File exists: iieout
/data/runs_geodyn/st/results/msis2/msis2_acceloff/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
/data/runs_geodyn/st/results/msis2/msis2_acceloff/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
/data/runs_geodyn/st/results/msis2/msis2_acceloff/DENSITY/st030914_2wk.goco05s
Loading data...
Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded
Compare the Coefficient of Drag¶
Details: - the \(C_d\) is an adjusted parameter in GEODYN. Each iteration, it is estimated and adjusted to improve the OD simulation - \(C_d\) is treated as a constant in GEODYN (for any iteration). The factor \(C_d\) varies slightly with satellite shape and atmospheric composition. However, for any geodetically useful satellite, it may be treated as a satellite dependent constant.
$ \bar`{A}_D = -:nbsphinx-math:frac{1}{2}` C_d \frac{A_s}{m_s} \rho_d v_r :nbsphinx-math:`bar`{r}_r $
\(A_D\) acceleration due to atmospheric drag force
\(C_d\) is the satellite drag coefficient
\(A_s\) is the cross-sectional area of the satellite
\(m_s\) is the mass of the satellite
\(\rho_d\)is the-density of the atmosphere at the satellite position
\(\bar{v}_r\) is the velocity vector of the satellite relative to the atmosphere, and
\(v_r\) is the magnitude of the velocity vector, vr .
Question: - What does a higher \(C_d\) actually indicate in a satellite drag scenario?
[7]:
labels = list(AdjustedParams_m1[5][SAT_ID]['0CD'].keys())
val_list = []
for i in AdjustedParams_m1[5][SAT_ID]['0CD'].keys():
val_list.append(AdjustedParams_m1[5][SAT_ID]['0CD'][i]['CURRENT_VALUE'])
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
from plotly.subplots import make_subplots
last_iter = 5
which_stat = 'CURRENT_VALUE' # 2 is the current val
labels = list(AdjustedParams_m1[5][SAT_ID]['0CD'].keys())
fig = make_subplots(
rows=1, cols=1,
subplot_titles=("Time Dependent Drag Coefficient (last iter)",
))
val_list = []
for i in AdjustedParams_m1[last_iter][SAT_ID]['0CD'].keys():
val_list.append(AdjustedParams_m1[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list,
name= m1,
mode='markers',
marker=dict(
size=10,),
), row=1, col=1,
)
val_list = []
for i in AdjustedParams_m2[last_iter][SAT_ID]['0CD'].keys():
val_list.append(AdjustedParams_m2[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list,
name= m2,
mode='markers',
marker=dict(
size=8,),
), row=1, col=1,
)
val_list = []
for i in AdjustedParams_m3[last_iter][SAT_ID]['0CD'].keys():
val_list.append(AdjustedParams_m3[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list,
name= m3,
mode='markers',
marker=dict(
size=8,),
), row=1, col=1,
)
fig.update_yaxes( title=r"$C_D (Unitless)$",exponentformat= 'power',row=1, col=1)
fig.update_xaxes( title="Date", row=1, col=1)
iplot(fig)
Orbit Sampled Density from Models¶
Details: - the density is sampled from the models along the orbit of the Starlette satellite
[8]:
data_nums = 1000
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
fig = go.Figure(data=[go.Scatter(x=Density_m1['Date'][:data_nums],
y=Density_m1['rho (kg/m**3)'][:data_nums],
name=m1,
mode='markers',
marker=dict(
size=4,
),
),
],
)
fig.add_trace(go.Scatter(x=Density_m2['Date'][:data_nums],
y=Density_m2['rho (kg/m**3)'][:data_nums],
name= m2,
mode='markers',
marker=dict(
size=2,),
),
)
fig.add_trace(go.Scatter(x=Density_m3['Date'][:data_nums],
y=Density_m3['rho (kg/m**3)'][:data_nums],
name= m3,
mode='markers',
marker=dict(
size=2,),
),
)
fig.update_layout(
title="Density along Starlette Orbit (for 4 minutes of first day)",
yaxis_title=r'$\frac{kg}{m^3}$',
xaxis_title="Date",
)
fig.update_yaxes(type="log", exponentformat= 'power',)
fig.update_layout(legend= {'itemsizing': 'constant'})
iplot(fig)
[9]:
data_nums = 100
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
fig = go.Figure(data=[go.Scatter(x=Density_m1['Date'][::data_nums],
y=Density_m1['rho (kg/m**3)'][::data_nums],
name = m1,
mode='markers',
marker=dict(
size=3,
),
),
],
)
fig.add_trace(go.Scatter(x=Density_m2['Date'][::data_nums],
y=Density_m2['rho (kg/m**3)'][::data_nums],
name= m2,
mode='markers',
marker=dict(
size=2,),
),
)
fig.add_trace(go.Scatter(x=Density_m3['Date'][::data_nums],
y=Density_m3['rho (kg/m**3)'][::data_nums],
name= m3,
mode='markers',
marker=dict(
size=2,),
),
)
fig.update_layout(
title="Density along Starlette Orbit (Every 100th datapoint)",
yaxis_title=r'$\frac{kg}{m^3}$',
xaxis_title="Date",
)
fig.update_yaxes(type="log", exponentformat= 'power',)
fig.update_layout(legend= {'itemsizing': 'constant'})
iplot(fig)
Residuals (obervation residuals)¶
Details: - in general, residuals are (observed minus computed) - obervation is from the observation file supplied by SLR tracking data - computed is what GEODYN is computing the location of the satellite should be -
these observation residuals show the WHAT EVEN ARE THE RESIDUALS
[10]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
import pandas as pd
fig = go.Figure(data=[go.Scatter(x=(Resids_m1['Date']),
y=Resids_m1['Residual'].values.astype(float),
name =m1 ,
mode='markers',
marker=dict(
size=3,
),
),
],
)
fig.add_trace(go.Scatter(x=(Resids_m2['Date']),
y=Resids_m2['Residual'].values.astype(float),
name= m2,
mode='markers',
marker=dict(
size=3,),
),
)
fig.add_trace(go.Scatter(x=(Resids_m3['Date']),
y=Resids_m3['Residual'].values.astype(float),
name= m3,
mode='markers',
marker=dict(
size=3,),
),
)
fig.update_layout(
title="Observation Residuals from Final Iteration",
yaxis_title="Residuals",
xaxis_title="Date",
)
fig.update_layout(legend= {'itemsizing': 'constant'})
iplot(fig)
Percent Difference in Density Output¶
[11]:
data_nums = 3
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
data_nums = 100
A = Density_m1['rho (kg/m**3)'][::data_nums]
B = Density_m2['rho (kg/m**3)'][::data_nums]
C = Density_m3['rho (kg/m**3)'][::data_nums]
diff1 = ((A-B)/B)*100
diff2 = ((A-C)/C)*100
diff3 = ((B-C)/C)*100
fig = go.Figure(data=[go.Scatter(x=Density_m1['Date'][::data_nums],
y=diff1,
name = m1+' and '+m2+' % diff',
mode='markers',
marker=dict(
size=3,
),
),
],
)
fig.add_trace(go.Scatter(x=Density_m1['Date'][::data_nums],
y=diff2,
name = m1+' and '+m3+' % diff',
mode='markers',
marker=dict(
size=3,
),
),
)
fig.add_trace(go.Scatter(x=Density_m3['Date'][::data_nums],
y=diff3,
name = m2+' and '+m3+' % diff',
mode='markers',
marker=dict(
size=3,
),
),
)
fig.update_layout(
title=r"Percent Difference in MSIS models along orbit of Starlette <br> Every 100th point",
yaxis_title='%',
xaxis_title="Date",
)
fig.update_layout(legend= {'itemsizing': 'constant'})
iplot(fig)
[12]:
data_nums = 5000
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
fig = go.Figure(data=[go.Scatter(x=Density_m1['Date'][:data_nums],
y=Density_m1['Z'][:data_nums],
name = m1,
mode='markers',
marker=dict(
size=3,
),
),
],
)
fig.add_trace(go.Scatter(x=Density_m2['Date'][:data_nums],
y=Density_m2['Z'][:data_nums],
name= m2,
mode='markers',
marker=dict(
size=3,),
),
)
fig.add_trace(go.Scatter(x=Density_m3['Date'][:data_nums],
y=Density_m3['Z'][:data_nums],
name= m3,
mode='markers',
marker=dict(
size=3,),
),
)
fig.update_layout(
# title="Density along Starlette Orbit (Every 15th datapoint)",
yaxis_title=r'$\frac{kg}{m^3}$',
xaxis_title="Date",
)
fig.update_yaxes( exponentformat= 'power',)
fig.update_layout(legend= {'itemsizing': 'constant'})
iplot(fig)
[13]:
# Trajectory_m2[7501001].columns
[14]:
# data_nums = 5000
# import plotly.graph_objects as go
# import numpy as np
# from plotly.offline import plot, iplot
# %matplotlib inline
# fig = go.Figure(data=[go.Scatter(x=Trajectory_m1[7501001]['Date'][:data_nums],
# y=Trajectory_m1[7501001]['HEIGHT'][:data_nums],
# name = m1,
# mode='markers',
# marker=dict(
# size=3,
# ),
# ),
# ],
# )
# fig.add_trace(go.Scatter(x=Trajectory_m2[7501001]['Date'][:data_nums],
# y=Trajectory_m2[7501001]['HEIGHT'][:data_nums],
# name= m2,
# mode='markers',
# marker=dict(
# size=3,),
# ),
# )
# fig.add_trace(go.Scatter(x=Trajectory_m3[7501001]['Date'][:data_nums],
# y=Trajectory_m3[7501001]['HEIGHT'][:data_nums],
# name= m3,
# mode='markers',
# marker=dict(
# size=3,),
# ),
# )
# fig.update_layout(
# # title="Density along Starlette Orbit (Every 15th datapoint)",
# # yaxis_title=r'$\frac{kg}{m^3}$',
# xaxis_title="Date",
# )
# fig.update_yaxes( exponentformat= 'power',)
# fig.update_layout(legend= {'itemsizing': 'constant'})
# iplot(fig)
[ ]:
[ ]: